Here we model microtubules with a finite width in order to properly account for the exclude volume interaction in the dense phases of cortical microtubule array. With high growth rate, a line segment model with zero width may lead to ultra-high density. The effective width of microtubule thus plays an important role to regulate microtubule density in such phases. A schematic illustration of the model is shown in Fig. 1. The basic building blocks of the model are rectangles with varying length l_j,k and width W, where the subscripts denote the k^th segment on j^th microtubule in the system. The very ends of micotubules are covered with semicircles with radius R = W/2. The interaction between microtubules is the hard rod potential, and any overlap of microtubules is strictly forbidden in the model. The growth of a microtubule is blocked if there is no non-overlapping space near the plus end. The nucleation of new microtubule can happen only in place where free-space for new microtubule with segment length Δl plus the circular caps is available.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. Schematic representation of an interacting kinetic growth model for microtubules with finite width W. The size of an individual microtubule j is described by a constant width W, and segment lengths lj,k, where the subscripts denote the kth segment on jth microtubule in the system. Every microtubule is covered with semicircle at both ends. The dash-dotted line indicates the central axis of microtubule, and the red arrow indicates the plus-end direction. In the cartoon, CMT2 will collide with CMT1 by further polymerization. αj is the colliding angle of the jth microtubule with its neighbor. The result depends on the choice of an critical zippering angle α0. For steep angle with αj > α0, the growth of microtubule is blocked and may trigger induced catastrophe; for shallow angle αj < α0, such as CMT3, it will zipper and grow along its colliding neighbor’s body axis. CMT4 represents a newly nucleated one in the free space.

Following previous studies, at present coarse-grained level, the single microtubule dynamics are sufficiently described by the growth/catastrophe/shrinkage/rescue on the plus end of the microtubule and shrinkage on the minus end.^[25,46,47] The growth speed k_g on the plus end is usually larger than the shrink speed k_sm on the minus end, otherwise the microtubule cannot exist for long time. The catastrophe rates on the plus end can be distinguished by spontaneous catastrophe rate k_sc and collision induced catastrophe rate k_ic. In the dilute limit, since collisions among microtubules are rare, we can expect that spontaneous catastrophe is more dominant than induced catastrophe.^[48,49] But here, since we are interested in the dense phase, collision happens so frequently that the induced catastrophe may dominate the catastrophe events.^[50] The catastrophe events will lead to the shrinkage of microtubule k_s on the plus end. Microtubules can resume growth with a rescue rate k_r. The typical values of the kinetic parameters used in the simulations are listed in Table 1.

Table 1.

Table 1.

Table 1. List of parameters and values in the simulation. . Parameters Simulated values kg plus-end polymerization rate 40 – 80 nm/s ks plus-end depolymerization rate 160 nm/s ksm minus-end depolymerization rate 40 nm/s kr rescue rate 0.01/s ksc spontaneous catastrophe rate 0.00–0.02 s−1 kic induced catastrophe rate 0.00–0.14 s−1 kn nucleation rate rate 3 μm−2⋅s−1 α0 critical zippering angle 0°–90° αc critical angle for collision induced catastrophe ∼40° Δl subunit length 10 nm L system size 8 μm T total simulation time ∼ 106 s W effective width 0–50 nm

Table 1. List of parameters and values in the simulation. .

During the growth, the microtubules may encounter each other. The growth of blocked microtubule will be stopped since no free space is available for the addition of a new segment. Microtubules are quite rigid, with typical persistence length of several millimeters.^[51] The small flexibility still allows it to deviate its original growth direction. In experiment, the so called zippering events are observed when two microtubules encounter, or even crossover if microtubule growth into the third dimension out of the plane.^[11,50] For simplicity, here we restrict ourselves to strictly 2D, thus we do not consider the crossover events. For zippering, it has been observed that such event depends on the relative angle α between colliding microtubules.^[11,52] We can define a critical α₀ to approximate the angular dependence. If the collision angle α>α₀, zippering is not possible. If the system allows induced catastrophe, then such large angle encounter would trigger such event with a rate k_ic. For α < α₀, zippering will take place, and the colliding microtubule changes its growth orientation following the direction dictated by the encountered microtubule (see Fig. 1).

We apply the event driven Gillespie algorithm to simulate our mechanical–chemical coupled dynamic system.^[29–31,53] In each simulation step, the time is forwarded according to , where r ∈ (0,1] is a uniform random number and a₀ = k_nL² + (k_gN₊ + k_sN_– + k_smN)/Δl + k_rN_– + k_scN₊ + k_icN_b is the sum of rates of all possible events. Here, N₊ and N_– are the numbers of microtubules in the growth-state and shrink-state on the plus-end, respectively, and Δl is the length of a discrete segment we add onto the microtubule. N_b is the number of microtubules which are blocked by encountering microtubules with the relative angle α > α_c, and the corresponding term is responsible for the induced catastrophe. α_c is the critical angle for collision induced catastrophe. Notice here that the term responsible for spontaneous catastrophe is k_scN₊, while for induced catastrophe is k_icN_b. N_b is the number of blocked configurations for microtubules, and in all simulation we have N_b ≤ N₊. In each simulation step, the algorithm will choose the stochastic event according to their relative weight. Once the algorithm decides a specific type of events (growth, shrinkage, nucleations, catastrophe, or rescue), it then randomly selects a specific microtubule with equal probability to realize such event. For growth and nucleation events, if we find the added segment with Δl overlap with the existing microtubule, we skip the growth/nucleation, keep the original configuration, and let the time increase. This leaves a gap at the tip of the blocked microtubule, which is a natural outcome of discretized nature of the segment.

The order of the system can be measured according the magnitude of weighted nematic order parameter tensor as defined in Refs. [15,30], which is given by

To see how the collective behavior of the dense microtubule array may be affected by the effective microtubule width and other dynamic processes, in the following sections we first use a simple tread milling model where the critical angle α₀ for zippering is set to 0°, and neglect catastrophe events on the plus-end of the microtubule. Next, we switch on the spontaneous and induced catastrophe events and demonstrate their relative importance in controlling the order–disorder transition. Finally, we vary the critical zippering angle α₀ and show the optimal regime to achieve nematic order.

To simulate crowded cortical microtubule array, we have to take into account of the effective width of single microtubule in the model. In vivo, the effective width of microtubule may differ from that of bare microtubule due to the complicated cortical environment and electrostatic interactions. First microtubules are strongly attached to the cell membrane, thus the membrane and associated proteins may mediate the effective interaction among microtubules. Moreover, the decoration of microtubule associated proteins may also change the effective width directly. In Fig. 2, we show how the effective width may change the order–disorder transition in the model system. In the phase diagram of Fig. 2, with the increase of the effective width W, the second nematic state emerges. Following previous studies, we denote the highly ordered nematic state with very long microtubules as N_I and relatively lowly ordered nematic state with short microtubules as N_II (see Figs. 3(a) and 3(b)). These properties have been well captured by the angular and length probability density function (PDF) of microtubules. As shown in Fig. 3(a), the angular distribution of the N_I phase is in a narrow range of angle around the super sharp peak, while the N_II phase although also has a dominant angle, indicating the emergence of nematic order, its distribution is broad. The length distribution also shows the singularity behavior of the N_I phase. The system shows a very broad length distribution until several thousand microns. The disorder phase and N_II phase, on the other hand, both show clear exponential decay of microtubule length.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. Phase diagram by changing (kg, W) parameters. Three phases are indicated, disorder isotropic phase (I) (pink region), nematic I phase with highly aligned and very long microtubules (NI) (blue region), and nematic II phase with relative short aligned microtubules (NII) (yellow region). Solid dots are samples for which simulations are performed. The star symbols represent the typical cases used in the analysis in Fig. 3, for I phase (W = 28 nm, kg = 45 nm/s) (see ESI-Movie1.mov), NI phase (W = 10 nm, kg = 65 nm/s) (see ESI-Movie2.mov), and NII phase (W = 40 nm, kg = 70 nm/s) (see ESI-Movie3.mov), respectively. Here the critical zippering angle α0 = 0. The spontaneous catastrophe and collision induced catastrophe are not taken into account in this phase diagram. Other parameters are given in Table 1. For every data point we performed one single thread simulation with time 105–5× 105 s.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Microtubule angular and length PDF for three typical cases as indicated by the star symbols in Fig. 2. (a) Angular PDF counted by subunits for three typical cases in three different phases. Isotropic phase I in green, sharply aligned NI phase in red, and aligned microtubules with broad angular distribution NII phase in blue. The distribution is obtained by time average of samples. Each sample is rotated with its dominant distribution at 90°. The vertical axis breaks at 0.1–0.9 to show the sharp red peak clearer (NI phase) which is about 20 times higher than the blue peak (NII phase). Inset: the vertical axis in log scale. (b) Length probability PDF of number of microtubules. Both isotropic (green) and NII (blue) phases demonstrate exponential decay in the length distributions. The NI phase shows a very broad length distribution, which indicates that the system is dominated by extremely long microtubules. The analysis is performed in the steady state for a period 2× 105 s with an interval of 103s. Inset: microtubule length distribution in NI phase with much broader range of length.

Although the symmetry of the N_II and N_I phases is the same, the transition between these two phases is usually discontinuous (first order-like), thus they could be well distinguished by studying how order parameter S varies with the control parameters. Here, we choose to fix the effective width W = 0 nm and W = 25 nm and increase the polymerization rate k_g so that we can expect a direct transition from isotropic phase I to N_I phase at W = 0 nm, and an intermediate transition into N_II phase at W = 25 nm. We measure order parameter S as defined in Eq. (1). As shown in Fig. 4(a), for W = 0 nm, the transition between isotropic phase and N_I phase is discontinuous, characterized by a coexistence regime around k_g = 44 nm/s. The discontinuity is a result of mutual enhancement of nematic order and increasing microtubule length. Without effective width, the transition can happen at relative low growth rate, with a direct jump between disorder and ordered states. For W = 25 nm, as shown in Fig. 4(b) the effect of crowdedness of the system is shown to change the transition from disorder to N_I phase. Although there still exist the mutual enhancement microtubule length and nematic order, the growth of microtubule is much easily blocked by other microtubules. With the increase of W, the microtubule occupied area will dramatically block the growth. Such strong steric effects thus delay the order–disorder transition, and can constrain the length distribution to be exponential even when the system becomes nematic ordered. As a result, the system can continuously go into the nematic N_II phase until very large polymerization rate, as shown in Fig. 4(b). When the order parameter S gets larger, the system finally breaks the blockage and dramatically elongates the microtubules with the formation of N_I phase. The large coexistence regime of N_II and N_I phases near k_g = 60 nm/s clearly indicates the discontinuity of such process.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. Isotropic nematic transitions for different effective microtubule widths. (a) Almost direct transition from disorder I phase to nematic NI phase for W = 0 nm with the increase of growth rate kg. The transition is discontinuous characterized by the coexistence regime near kg = 44 nm/s. (b) For W = 25 nm, the transition from disorder I phase to NII phase, then jumping into NI phase with the increase of growth rate kg. The NII and NI phase transition is discontinuous, with a bistable regime near kg = 60 nm/s. The order parameter S is obtained by time and sample average (20 samples). Errors bar for 20 independent samples are given and most are within the symbol. Red circles represent that simulation samples undergo phase transition process and do not reach the steady state within the simulation time ∼ 5× 105 s, giving the strong discontinuous nature of the transition, which is not included for averaging. Other parameters not listed here are the same as those in Fig. 2.

As a result, in the phase diagram of Fig. 2, the border between N_I and N_II/I phases is first order-like (discontinuous), and the border between I and N_II phases is second order-like (continuous). The interception point of these two transition lines is the so called tricritical point. A similar topology of the phase diagram also exists in a model without the inclusion of the effective microtubule width.^[15] In such case, the N_II phase will emerge if the system has a very large nucleation rate to produce lots of short microtubules in the model with W = 0 nm. Here, we have kept the nucleation rates unchanged, even effectively reduced with an effective width, but the inclusion of effect width increases the crowdedness of the system without further sharpening of the order. This plays the similar role as the nucleation of short disorder microtubules in previous model.

With the inclusion of the effective width, we now check how the catastrophe events will be coupled to the finite width microtubule system. The collision induced catastrophe depends on the blocked configurations among growing microtubules. In the dilute limit of the microtubule system, with negligible collisions, the spontaneous catastrophe plays much important roles. It regulates the length distribution and thus also the order property of the system. In the dense phase, the induced catastrophe becomes much important due to the frequent collisions among microtubules. We are interested here to see how different catastrophes may effect the self-organization of the dense phase.

For this purpose, we choose W = 25 nm (approximately the bare width of microtubule) and k_g = 70 nm/s with rescue rate k_r = 0.01 s^–1 (see Table 1), where the system stays in N_I phase without introducing catastrophe. To compare the difference of two catastrophe mechanisms, we focus on two cases where spontaneous catastrophe and induced catastrophe are introduced separately. As shown in Fig. 5(a), both catastrophe events can drive the system decay from N_I to N_II and then disorder. The catastrophes trigger the depolymerization process k_s on the plus-end, and the microtubule thus will shrink on both ends, and the length is shortened efficiently. This will reduce the microtubule length and lead to less ordered state. The transition from N_I to N_II is still discontinuous (first order-like), and N_II to isotropic is continuous (second order-like). The transition from N_I phase to N_II is more sensitive to spontaneous catastrophe k_sc. In N_I phase, the system is extremely crowded, however, the microtubules are highly aligned. Perfect alignment dramatically reduces the collisions among microtubules, which makes the transition less sensitive to induced catastrophe rate k_ic.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. Spontaneous catastrophe and induced catastrophe in regulation of order–disorder transition. (a) The transition from NI phase to less ordered NII phase and disorder with the increase of spontaneous catastrophe rate ksc (red square) or kic (blue dot). (b) Same data as (a) but the horizontal-axis is mapped to the time averaged frequency of real catastrophe events (RCE) for all the samples. The critical angle for collision induced dynamic catastrophe is αc = 40°, the growth speed kg = 70 nm/s, and effective microtubule width W = 25 nm. The simulations start from NI state as obtained in Fig. 4(b) at kg = 70 nm/s. Each point is averaged with 5 samples with time interval 1 s, and average simulation time longer than 105s.

Since k_ic needs the collision process to trigger the catastrophe, its value does not represent the frequency of catastrophe events that actually take place in the system. We can further compare the effects of both catastrophes by the frequency of real catastrophe events (RCEs). We count the average frequency of RCEs in simulations for both cases in steady state. The statistics are shown in Fig. 5(b).

There are two important features in such representation. First, compared with Fig. 5(a), the significant overlap of N_II and N_I phases for induced catastrophe disappears. It signifies that the N_II and N_I phases have very different RCEs even in the coexistence regime. For given k_ic in the coexistence regime in Fig. 5(a), the RCEs in highly ordered N_I phase are almost 0, while for relatively lowly ordered N_II phase, there exist quite significant collisions among microtubules, thus quite large RCEs are obtained.

Second, for given RCEs, the induced catastrophe still demonstrates significantly higher nematic order than the case for spontaneous catastrophe. The spontaneous catastrophe events only depend on the number of microtubules that are in the growth state, as we can see from the term k_scN₊ in Eq. (1). However, the induce catastrophe events further require that these microtubules are blocked by neighbors (k_icN_b in Eq. (1)). Since micortubules that are not aligned with their neighbors are much easier to be blocked during the growth, collision induced catastrophes are much efficient to recognize and clean the discordant microtubules than spontaneous catastrophe, thus allowing larger nematic order at the same RCEs.

Zippering also depends on the collision among microtubules. In an active dense phase, we expect it would play an important role in controlling the order–disorder transition. One important parameter controling the zippering events is the critical angle α₀. When the colliding angle α is less than α₀, two microtubules will zipper. However, when colliding angle α is larger than α₀, we expect it is going to block the growth, and even can trigger collision induced catastrophe or crossover events. As we show above, catastrophe events basically are going to reduce the order in the system monotonously. Its effects can be compensated by increasing the growth speed. Here to investigate the effect of zippering, we set catastrophe to 0 s^–1. We also restrict ourselves strictly in 2D, thus do not consider the crossover events.

Zippering is considered a promotion effect of order in the system.^[11] It helps to change dis-aligned microtubule groth along with their neighbor thus effectively introduces an alignment interaction. On the other hand, if zippering can happen very easily for large α, the microtubules essentially become very flexible. For flexible polymer, nematic order will disappear. To compromise the easiness of alignment and rigidity, we may expect an optimal critical zippering angle α₀ for nematic alignment.

In Fig. 6(a), we scan the (α₀, k_g) parameter space. It is clear that for large critical angle α₀ > 70°, the system becomes disorder even when k_g is large. For small α₀, at moderate growth rate 46 nm/s < k_g < 55 nm/s, the system is also in a disordered state. The N_II phase exists only for an optimal α₀ region. Large α₀ drives the system into disorder. For small α₀, the N_I phase becomes dominant when k_g is large. In Fig. 6(b), these properties are further shown by an explicit plot of order parameter S versus α₀. In Fig. 6(a), we also notice that, near k_g = 55 nm/s, the reentrant behavior between N_I and N_II phases indicates that the zippering can modify the ordering in a non-trivial manner. In such condition for α₀ = 0°, the microtubules are always blocked during collision, and the system is locked in N_II phase. However, when a small α₀ is introduced, it allows the colliding microtubule to grow if the colliding angle α < α₀, and this suddenly drives the system into highly ordered N_I phase. Further increase of α₀ however makes the microtubules too flexible to maintain the N_I phase.

Fig. 6.

Figure Option ViewDownloadNew Window

Fig. 6. The effects of critical zippering angle on growing microtubule phase behavior. (a) Phase diagram with the variations of critical zippering angle α0 and kg. The phase diagram is colored the same way as in Fig. 2. (b) The variations of time averaged order parameter S with the change of α0 and kg. The effective microtubule width used here is W = 25 nm, and kic = ksc = 0 s–1. The average is performed in the steady state with a sampling interval 1 s and total average period 105–5× 105 s.

The optimal zippering angle indicates that zippering does not simply enhance the alignment. Although it enables the alignment of the growth direction of microtubules, it leaves those dis-aligned segments unchanged. Such dis-aligned segments can only be eliminated by depolymerization on the minus end. The catastrophe events on the plus-end would not help since it starts from the plus end and will eliminate the aligned part first. For large α₀, the zippering events happen irrespective of the colliding angle, which leads to frequent large angle bending of microtubules and drives the system into complete disorder.

In our simulation, by systematically varying the critical zippering angle α₀, we find optimal zippering angle around 30° to 40° for N_II phase. This optimal angle is consistent with the experimental determined zippering value. In vivo, various conditions may help regulating the critical angle for zippering. Since microtubules have large persistence length, the most important factor in control might be the microtubule associated proteins that fastens the microtubule on the cell cortex.^[36,44,54] Higher density of these binding proteins might allow small critical zippering angle. Moreover, the plus-end microtubule associated proteins might also play a role in the regulation of zippering. It would be interesting to investigate if evolution has allowed plant cell to exploit such possibility to do optimization for achieving ordered cortical microtubule array.